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ABSTRACT 


Plate problems often involve dynamic disturbances by time dependent external 
forces or displacements. The purpose of this study is to investigate the effect of dynamic 
loading on propagation of cracks in a plate. As it is not possible to obtain closed form 
solutions for dynamic fracture mechanics problems of practical interest, one has to make 
use of numerical method. 

The existing force release models to simulate crack propagation have certain 
drawbacks giving raise to large oscillations in the solution. In the present work a new 
model is developed to study high speed propagation of cracks. In this model instead of 
reducing the crack tip node force, an additional one-dimensional elastic element is 
attached to the crack node, and the stiffness value is gradually reduced as the crack 
proceeds to next node. 

The FE code is first validated with a static problem and after validation it was 
used to solve quasistatic and dynamic crack propagation problems and then extended to 
solve the inverse problem. The energy release rate due to the crack propagation with 
respect to time is calculated for different propagation speeds and found to be more 
consistent with practical results. 
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CHAPTERl 


INTRODUCTION 

1.1 E^TRODUCTION: 

The dynamic fracture phenomenon can be characterized by various dynamic states 
of crack tip. The dynamic states of crack tip are induced by impact loading applied to the 
cracked solids, or by fast motions of the crack tip itself. Dynamic loads may be created by 
moving vehicles, wind, gusts, shocks or blasts, unbalanced machines, wave impacts, seismic 
disturbances etc. For solving these problems often the time variation is neglected and they 
are treated as static or quasistatic problems. Structural members which are subjected to 
cyclic loading such as vibrations, microscopic imperfections in the material give rise to 
microscopic cracks which in turn grow into cracks of significant size over passage of time. 
At this stage further growth of cracks stably or unstably leading to the component failure 
depends not only on the material characterising but also on the nature of loading, i.e static 
or dynamic. For many dynamic problems it is impossible to find closed form solutions, 
especially if they involve the aspect of fracture. A number of analytical solutions to 
problems of dynamic fracture shedding important light on the basic phenomenon have 
been obtained in the past two decades. These solutions are however limited to simple 
cases of loadings and of unbounded plane bodies. Usually, the interactions of the stress 
waves emanating from the crack tip and or those reflected from the boundaries of a finite 
body make the analytical solutions invalid. Thus, it often becomes mandatory to analyze 
dynamic crack propagation in finite solids using numerical methods. 
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The subject of the mechanics of dynamic fracture can be broadly described as that 
branch of mechanics of solids, where in the effect of material inertia and stress waves on 
crack play significant roles. One method of classifying dynamic fracture problems is as 
follows [Nishoka and Atulri, I986j; 

(1) Solids containing stationary cracks subjected to dynamic loading 

(2) Solids containing dynamically propagating cracks under quasi-static loading 

(3) Solids containing dynamically propagating cracks under dynamic loading 

A knowledge of time-dependent asymptotic stress and displacement fields near the 
crack tip, and their "strength "as quantified, for instance, by the "stress intensity factor (SIF) 
in linearly elasto-dynamics, is basic in understanding the process and nature of dynamic 
fracture in solids. A direct approach to laboratory evaluation of dynamic fracture 
toughness and crack arrest toughness, is based on the laboratory measurements of the 
instruments dynamic stress fields close to the propagating crack tip using photo elastic 
method. Such direct measurements, however, are generally difficult. To overcome these 
difficulties, and simplify the measurements necessary, hybrid experimental-numerical 
methods are often preferred. This is one reason that the advancement of the state of 
science of dynamic fracture mechanics relies heavily on simultaneous advances in 
computational methods. 

Dynamic fracture mechanics concerns (i) the onset of crack growth under dynamic 
loading (ii) dynamic crack propagation in stressed solids. The 

conventional linear elasto-plastic (or) quasistatic fracture mechanics applies only up to the 
end of stable growth under quasistatic loadings and assumes that the onset of unstable 
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crack propagation renders the structure useless. Yet, the prevention of crack growth 
initiation itself may be too conservative or too costly an objective for some structural 
designs, and also catastrophic failures caused by unstable crack growth are obviously 
intolerable. In these cases, the assurance of crack arrest, as a second line of defense, is 
essential. This concept has been investigated for design methodologies assuring the 
integrity of nuclear pressure vessels under thermal shock conditions, LNG ship hulls, and 
gas transmission pipelines. 

The fundamental frame work of the subject of dynamic fracture mechanics relies 
primarily on solutions for dynamic behaviour of solid containing cracks these solutions are 
characterized by (i) stress wave integrations (ii) the need to account for kinetic energy in 
the global energy balance of the fracturing body and (iii) inertia effects of the material. 
The inertia effect may arise due to two reasons (i) rapidly applied load on a cracked solid 
and (ii) rapid crack propagation. In the first case the influence of the load is transferred 
to the crack by stress waves through the material. In the second case the material particles 
on the two crack faces displace with respect to each other, as the crack advances. The 
inertia effect, in the first case, is considered significant when the time taken to load the 
specimen to maximum value is small as compared to the time required for a characteristic 
dimension of the body. In the second case, the inertia effect should be accounted far 
whenever the crack velocity is a significant fraction of the characteristic velocity. 

) 

In the past two decades or so, a number of analytical solutions, which 
provides a useful understanding of dynamic crack behaviour, have been obtained. These 
analytical solutions are however, limited to cases of simple loadings and unbounded plane 
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bodies. Moreover the stress wave introductions, which play an important role in dynamic 
fracture mechanics. Usually render the analytical solutions intractable. Therefore, the use 
of numerical methods is often indispensable for the analysis of cracks in finite solids. 

Numerical methods in fracture mechanics were initially appraised in 1978 at that 
time, a comparison of the finite element and finite difference methods to the following 
conclusions. The finite element method was more suitable, for the analysis of stationary 
cracks under dynamic loadings. Due to the fact that the relevant singularities can be 
modeled in the crack tip elements. In large measure, this is the result of the development 
of finite element procedures to model propagating singularities and of path-independent 
integrals which characterize the strength of the fields near propagating crack tip. 


1.2 LITERATURESURVEY: 


Though considerable amount of work have been done to study the fracture 
phenomena through numerical methods, only large rectangular double cantilever beam 
specimens were given special interest Owen and Shantaram (1977) started the case of finite 
element method to study dynamic crack growth. They studied double cantilever beam 
specimen and pipeline problems under transient loading. The crack was advanced from 
one node to the next node, when the stress at the gauss point nearest to the crack tip 
exceeded a certain value and then damping coefficient was made zero at the released node. 
The crack propagation history was simulated for the given loading conditions. 
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Nishoka and Atluri (1982a) investigated the crack propagation and arrest in a high 
strength steel DCB specimen using moving singular dynamic finite element procedure. An 
edge crack in a rectangular DCB specimen was propagated by inserting a wedge. The 
results were compared with experimental data of caustics. In another work, Nishoka and 
Atluri (1982b) presented the results of generation and prediction studies of dynamic crack 
propagation in plane stress and plane strain cases. The studies were conducted by using 
FEM, taking into account stress singularity near the crack tip. The variation of dynamic 
stress intensity factor with time and the variation of dynamic fracture toughness with 
velocity were studied and compared with available experimental results. 

Nishoka and Atluri (1986) gave elaborate information on analysis of dynamic 
fracture using FEM. The most common ways to deal with the crack tip region are to 
simulate crack growth through gradual release of element nodal forces or embedding a 
moving element in which the interpolation functions are determined by the continuum near 
tip fields at the crack tip in the mesh, and J-lntegral considerations. Following spacial 
description the differential equations in time for the node point variables must be 
integrated. Because the dynamic fields associated with rapid crack growth are rich in high 
frequency content, small time steps are required for accuracy. Because of natural 
restriction to small time steps, the authors report that many times combination of 
conditionally stable explicit time integration scheme and diognalized mass matrix are found 
to be accurate. 

Cronech and Williams (1987) used dynamic and generation mode finite element 
program to analyze different specimen geometries. Chiang (1990) presented a numerical 
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procedure based on eigen functions to determine the dynamic stress intensity factor of 
crack moving at steady state under antiplane strain condition. An edge crack problem and 
a radial crack problem was solved using traction and displacement boundary conditions 
separately. It was shown that the dynamic effect is relatively insignificant for low crack 
propagation speeds provided that specimen size fairly is large. 

A detailed analysis of the various time integration schemes from the point of view 
of accuracy and speed, were reported by (Scran et al., 1990; Zienkiewicz and Taylor, 1991, 
Bathe, 1990). It was concluded that, even though central differences scheme is the best in 
terms of computational speed and cost, Newmark’s constant average scheme gives better 
accuracy and offers unconditional stability, whereas, the time step size has negligible 
effects when central difference time integration scheme is employed (Bazant and Celep, 
1982, 1983). 

Thesken and Gudmundson (1991) worked with an elasto-dynamic moving element 
formulation incorporating a variable order singular element to enhance the local crack tip 
description. Kennedy and Kim (1993) incorporated micropolar elasticity theory into a plain 
strain finite element formulation to analyze the dynamic response of the crack. Material 
with strong micropolar properties were found to have significantly lower dynamic energy 
release rate than their classic material counter-parts, Wang and Williams (1994) 
investigated high speed crack growth in a thin double cantilever beam specimen using 
FEM. The cantilever end was loaded under an initial step displacement of 1 mm and then 
pulled with a constant velocity the crack was propagated at assumed speed by gradual node 
release technique. Large dynamic effects were observed because of wave reflections within 
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finite specimen size. The reflection and interaction of the waves from the free boundary 
were avoided by limiting the duration of study. 

1.3 PRESENT WORK: 

In the present work a new model is developed to study highspeed propagation of 
cracks in a 2-dimensional plate, for simplicity the crack is symmetric and is under model 
propagation. In this model instead of reducing cracktip node force, an additional one- 
dimensional elastic element is attached to the crack node , whose stiffness value is 
gradually reduced to zero as the crack proceeds to next node. This model is expected to 
avoid the oscillations in solution. The present work makes use of some parts of the code 
developed for the dynamic analysis of two dimensional problem by Kishore, Kumar and 
Verma [1995]. 

Chapter 2 deals with the finite element formulation and evaluation of fracture 
parameter and explains the details of crack propagation models. Chapter 3 discusses the 
new model developed in the present analysis. Chapter 4 deals with the validation of code 
and results for the test problems. Chapter 5 deals with conclusions of the present 
work and the scope of further work. 
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CHAFrER2 


FEM FORMULATION AND CRACK MODELLING 


2.1 FINITE ELEMENT FORMULATION: 


The hyperbolic partial differential equation governing elastic wave propagation in 
homogeneous,, anisotropic and linear elastic solid is. 


The total potential of a linear elastic body, neglecting body forces and using D’Alembert’s 
principle, is given by 


“ = L - 1 


uf. ds 


L 




dv 


( 2 . 2 ) 


where the first term is the strain energy, second term is the workdone by the surface 
traction tj and the last term is the work done by the inertia forces. The integrations are 
carried out over the volume of the body v or surface s, whichever is applicable. 

The finite element solution involves the discretization of the domain into suitable 
elements, approximation of the field values interior to the 

elements in terms of its nodal values through the shape functions of the chosen elements 
and the determination of the nodal values through the minimization of the above energy 
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functional. Invoking the stationarity of tt, namely, Stt = 0, the following equation may 
be obtained 


L 




- f SWjf, ds + f du.pu. dv = 0 

J S V V 


(2.3) 


Approximating the field values (u) within an element in terms of its nodal values (U) as 


u^(x,y,z,t) = N%x,y,z)U'^(t) 


(2.4) 


From this following strain-displacement relation is obtained 

e''(x,y,z,t) = B‘(x,y,z)U'^(t) (2-^) 


The constitutive equation is 

r^(x,y^,t) = D‘ e‘(x,y^,t) 


( 2 . 6 ) 


N® is the displacement interpolation matrix, B® is the strain-displacement matrix and D® 
is the material property matrix. 

Substitution of Eqs. (2.4) to (2.6) into (2.3) gives the equation of motion in the matrix 
form as. 


[M] 


c^U 

dt^ 


H- [jq [[/] = [FI 


(2.7) 


where [Mj and [K] are the system global mass, stiffness matrices are obtained by 
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assemblage of elemental stiffness, mass matrices and [FJ is external load vector. [Uj is 
global displacement vector of finite element assemblage. 

m = nMY, [K\ = E[Jq^ [F] = 


where [Mj®, [Kj® and [F]® are the elemental mass, stiffness and load matrices 
respectively. In the wave propagation analysis the body is usually descritized using 
simple elements such as 3 -noded triangular or 4 -noded quadrilateral element to avoid 
spurious reflections at element boundaries. In the present work 4 -noded element is 
used. Using the general terminology of the elemental matrices can be briefly explained 
as follows: 


[Ml' = f pim^lNldv 


(2.8) 


[K]^ = f 

Jv, 


(2.9) 


[F\^ = f minds 


( 2 . 10 ) 


where 

[NJ is set interpolation functions by which one can obtained displacements at any spatial 
point from nodal values. 

[BJ is set of derivative of interpolation functions. 

[D] is the material matrix. 

[t] is the traction vector. 
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The equation 2.7 can be solved either by time integration or by mode 
superposition of which the former is preferred in the wave propagation problem. In this 
integration scheme, there are many different methods, which can be classified as 
"explicit" or implicit" whose advantages and disadvantages of each method are given in 
detail in Bathe (1990). In the present work Newmark integration method for time 
variable is used. The details of the Newmark method are given in Bathe (1990). The 
choice of time step At for time integration is an important factor for an accurate solution. 
An optimum choice of time step is At = d/Cj in which d is the smallest element size, 
and Cj is dilatation wave velocity. 

2.2 EVALUATION OF FRACTURE PARAMETERS 

The crack effects are usually characterised by stress intensity factors (Kj, Ku, 
Kiii),energy release rate (Gj, Gjj, Gjjj) and J-Integral. While the SIF’s require 
special (moving) elements to be used, the other two give reasonably accurate results with 
ordinary meshes. These two parameters are described in the following sections. 

2.2.1 Energy Release Rate 

Griffith proposed that a crack in a body will not extend unless energy is released 
in the process to overcome the energy needs of forming two new surfaces, one below and 
one above the crack plane. As the crack in the plate advances, the plate becomes less 
stiff, consequently, for the case of the plate with ends held rigidly, the stress within the 
plate decreases, and the strain energy stored in the plate is reduced. The energy thus 
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released is available for the crack to grow. Most of the energy released, as the crack 
advances comes from those parts of the plate adjacent to the cracked planes. There are 
two important quantities involved (i) how much energy is released when a crack advances 
and (ii) minimum energy required for the crack to advances in forming two new surfaces. 
The first quantity is measured in terms of a parameter energy release rate, denoted by 
‘G’. Thus the energy release rate is defined as energy release per unit area extension 
of crack growth if there is a virtual crack growth, i.e, energy equal to G would be 
released from the system per unit extension of area. The energy requirement for a crack 
to grow by a unit area is called crack resistance. Crack resistance is sum of two energies 
required to (i) form new surfaces and (ii) cause anelastic deformation. Both the 
available energy release rate and crack resistance are important to study the possibility 
of crack becoming critical. Obviously, the available energy release rate must be greater 
than crack resistance for a crack to have a chance of grow. If the energy release rate 
exceeds the crack resistance, the crack starts to grow at high speeds. 

Mathematical Formulation: 

As the crack advances the following phenomena happen: 

1. Energy in the component decreases 

2. Stiffness of the component decreases 

3. Energy is consumed to crack two new surfaces. 

Formulation for energy release rate is carried out invoking the principle of conservation 
of energy. First, consider quasistatic crack growth case. Let the incremental increase 
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in the crack area be AA. This crack growth is achieved by an incremental external work, 
A done by the external forces and the strain energy within the body of the 

component increases by AU then available energy, GAA provides the energy balance as 
follows: 

GAA=^W^-^U (2.11) 

a . ( 2 . 12 ) 

oA 

G = (2.13) 

dA 

where tt is potential energy 

If the plate of uniform thickness B then AA can be expressed as B da. 
da is increment in crack length. 


In dynamic case, as crack moves rapidly some energy is being consumed to impart 
kinetic energy to cracked portion of the body and to generate stress waves. Therefore G 
in the dynamic case can be far different from quasistatic case. Thus Eq.3. 1 can be written 
as 

GAA = AW^-AU-AT (2.14) 


where AT is increment in kinetic energy in the body 


— ( 


dU 

dT 

b\ 

da 

da 

da 


(2.15) 
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in terms of FE analysis the above terms can be written as 


u = 

2 

T = 

where 

= external work 
U = strain energy 
T = kinetic energy 
[RJ = external load vector 
lU] = displacement vector 
[U] = velocity vector 
[K] = Stiffness matrix 
[Mj = mass matrix 

2.2.2 J-integral: 

Under approximate assumptions of material homogeneity, the strength of the 
crack tip field is governed by an integral evaluated over a path that is far removed from 
the crack tip. Since the stress and displacement data are evaluated away from the crack 
tip, these values are relatively insensitive to the linear details of modelling of crack tip 
region. J-integral is used for computing elasto-static fracture is mode I. Also the far 
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field integral can be calculated with reasonable accuracy using relatively course finite 
element models of the structure. 


Atluri (1982) and Nishioka and Atluri (1983) have under taken the studies of 
several path independent integrals in the analysis of crack growth initiation and 
propagation of cracks in elastic or inelastic materials under quasistatic and dynamic 
conditions. In the present work J is defined as follows 


J = 


lim 

c-0 


P«, “u dy 


(2.16) 


where W is strain energy density, nj is the direction of unit outward normal, t* is the 
surface traction and definition of the paths r^, r, r^. and the volume V, are shown 
in Fig. 2.1. The value of J is independent of the choice of r of Fig. 2.1 only under 
steady state crack growth conditions. The path is held stationary as the crack tip is 
extended in a self similar manner. 
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2.3 MODELLING OF CRACK PROPAGATION: 


To simulate the crack propagation in solids two different concepts of 
computational modelling have been in use, ie. the stationary mesh procedure and moving 
mesh procedure. Out of these two moving mesh procedure is the name implies, involves 
the changing of mesh for each iteration as the crack propagation. This is computationally 
very intensive. Stationary mesh procedure does not require change in the mesh in 
general, and requires only change in boundary and loading conditions. This procedure 
is explained in detail in the following paragraphs. 

Stationary mesh procedure: 

In the "stationary mesh" procedure of modelling linear elastodynamic crack 
propagation, when the node spacing is cAt (c being crack velocity), one only has to 
release the constraint on displacement at tlie old crack tip. If one also applies equal and 
opposite nodal forces to eliminate reaction forces at the old crack tip, the crack surface 
becomes over relaxed. 

As the crack propagation velocity is lower than the wave velocities, in general, 
the crack-tip occupies positions in between the nodes at various time step of the 
numerical integration. Thus if a simple i.e., crack propagation procedure by advancing 
tip from node to next node is used, then the crack-tip either stay at one node or jump 
from one node to next node in time At in the computer simulation. This causes the 
crack tip movement to be irregular or discontinuous, serious violations of the kinematics 
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and fraction condition ahead of the actual crack-tip will result. Also sudden increase of 
crack length by the release of constraints on displacement induces spurious high 
frequency oscillations in the finite element solutions. To overcome such difficulties, 
several algorithms have been suggested in literature to release the node gradually. 
Suppose that the actual tip is located at‘C’ is between the finite element nodes B and D 
as shown in Fig. 2.2. The lengths segments BC and BD are b and d respectively. The 
holding back force, F, at node b is gradually reduced to zero as the crack tip reaches to 
the node D. The various scheme available to decrease the force to zero as follows : 

(i) Malluck and King (1978) suggested the release rate based on constant 
stress intensity factor; that is 


JF^ 

Fo 



(2.17) 


Where F^ is original reaction force when the crack tip was located at node B. 

(ii) Pydham et.al (1978) suggested the release rate based on constant energy 
release rate. 



(2.18) 


(iii) Kobayasi et.al (1978) suggested the linear release rate based on no 
physical argument other than pure intuition. 


Fo 



(2.19) 
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For having gradual and smooth opening of crack Kishore, Kumar and Verma (1993) 
tried an alternative method. The holding back force at the crack tip B is linearly 
decreased to zero when the crack reaches the end of the next element. Thus, when the 
crack tip goes beyond node B then 

( 2 . 20 ) 

ly 



Where Fjjg is reaction at node B, when the node was closed; b is the crack extension 
and d, dj are the element lengths as shown in Fig. 2.2 and when crack moves beyond 
node D to point Dj 



d+b.\ 

1 1 

d+d, 

\ 1/ 


( 2 . 21 ) 




( 2 . 22 ) 


Where Fjjjj is the reaction at node D, when the nodes was closed; b, b 2 are the crack 
extensions and d, dj, and d 2 are the element lengths as shown in Fig. 2.2. 


In the present work instead of reducing the force a one-dimensional elastic 
element is attached to the crack node whose stiffness is reduced gradually as the tip 
advances to the next node. This model explained in detail in Chapters. 
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Fig. 2.2 Crack opening scheme in force release model 
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CHAPTERS 


STIFNESS RELEASE MODEL 

3.1 INTRODUCTION 


In the previous Chapter force release models to simulate crack propagation were 
described. These models when applied to high speed crack propagation have the 
following drawbacks; 

(1) The time increment At, is usually taken such that it takes 15 to 20 iterations for 
the cracktip to cross one element. The amount of force necessary to keep the 
nodes together is determined and a factor is used to proportionately decrease this 
force as the crack advances. Thus, at subsequent times, though the problem is 
highly dynamic, the current calculations were based on the force calculated to 
keep the nodes closed 15 to 20 time steps before. This will cause discripancy as 
the crack advances further from one element to the next one causing large 
oscillations in solution. 

(2) In the case of decreasing energy release rate and possible crack arrest in between 
the nodes, freezing the force at particular value in the dynamic problem will lead 
to a non acceptable solution. This aspect acquires importance in the case of 
inverse problems. 

In the present work to avoid the above difficulties a modified model has been 
developed. 
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(C) 


Fig. 3.1 Crack opening scheme in stiffness release model 
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3.2 STIFFNESS RELEASE MODEL 


In this model instead of applying a force at the opening node a one dimensional 
elastic element is added (Fig. 3.1) whose stiffness is reduced gradually as the crack tip 
advances to the next node. This additional stiffness lets the crack to be opened in a 
controlled manner. It is clear that the amount of additional stiffness required is very 
large (theoretically infinite) when the crack tip is at the begining of element and zero 
when the crack tip reaches the other end. An expression for additional stiffness to 
be added is derived with reference to quasistatic crack propagation. 

Fig. 3.1(a) shows a symmetric part of a deformed plate in which crack is upto the 
node B. Fig. 3.1(c) shows the configuration of the plate when the crack has completely 
propagated to the next node C. Fig. 3.1(b) shows the modelling of the plate with 
partially propagated crack (somewhere between B and C) by the addition of a one- 
dimensional linear spring element whose stiffness is Kg. Let Kq be the stiffness of the 
node B of the original plate without additional spring (i.e., force required to cause unit 
displacement of node B), Uq be the displacement of the node when the crack tip reaches 
the next node and the element completely opens. Then, Kg can be found from the 
equilibrium equation: 




(3.1) 



(3.2) 


Now to express as function of crack size, it is necessary to know ’u’as a function of 
crack size. For this purpose energy release rate is used. At any stage, the crack can be 
extended by a small amount by reducing Kg by an amount A Kg. This reduction in the 
spring stiffness causes a decrease in the total energy of the plate. As this is the only 
energy loss in the system, this can be directly used while evaluating energy release rate. 

The energy in the spring can be written as 


E = - Ku'^ = - 
2 * 2 “ 


u 


- 1 




(3.3) 


The energy release rate is written as, 

G = — (3.4) 

LA 


where A A is increment in crack area. 

Let a non-dimensionalised parameter for crack length, ’a’ be defined as, 
a = crack length (within the element) / d 
where d is the distance between two nodes. 

Therefore A A can be written as, 

AA = B d Aa 

where B is the thickness of plate. 

Then energy release rate, G becomes 
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(3.5) 


G = 


j. 

2 


^0 (-“o) 


1 du 
Bdda 


To determine the variation of u as a function of ’a’ it is necessary to know the variation 
of G. As it can be assumed that within the element G is linear function of crack length, 
let it be assumed that 


G = Cj + qa 


(3.6) 


Where Cj and Cj are constants to be determined later. 
Substituting G in the Eq. (3.5) 


+ CjO 


1 

2 Bd \ da j 


(3.7) 


(Cj + Cjfl) Bd da = du 


(3.8) 


Integrating on both sides 


Bd 


qa + + q 


-1 K^u^u 


(3.9) 


where Cq is integration constant 


9 I 

„ _ 2d qa + — — + q 

Vo I 2 


(3.10) 
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By applying initial condition 
using the other end condition 


u - 0 when a = 0, it can be seen that Cq = 0 

u= Uq when a = 1 in Eq. (3.10) it can be obtained. 


Un 


^ 0*^0 


Bd 



(3.11) 


From Eq.(3.10) and Eq. (3.11) it can be seen 

n 2C^a + 

Uq 2Cj + Cj 


Eq. (3.12) can be simplified as 

— = a + C (a-a^) 

Mn 


where 


ICy+C^^ 


Therefore, finally Eq. (3.2) can be written as, 


Ks = 





(3.12) 


(3.13) 


(3.14) 


where f = u / Uq given by Eq. (3.13) 
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As can be seen that Kq has to be determined for the given plate under quasistatic 
conditions and given boundary conditions. However, it may be inferred that the near 
boundaries affect the value of Kq than the far off boundaries. Suitable value of ’C’can 
be determined by trial and error to obtain continuous energy release rate in between 
elements. Usually it is a very small value. In the present 2-D problem it is of the order 
of 0.05. 

3.3 PREDICTION OF CRACK INITIATION, PROPAGATION AND ARREST: 

The present method can be easily applied to predict the crack propagation history 
under the known loads, if crack propagation resistance of the material is known. That 
means that fracture resistance of the material, G as a function crack propagation velocity 
is supposed to be a known input data (Fig. 3.2). At each time step the avaliable energy 
release rate is determined and it is compared with the material fracture resistant data 
to predict the crack propagation velocity. This is used to determined the crack 
propagation length in the th6 given time interval and propagated accordingly. This 
procedure continues till crack arrest. 
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CHAPTER 4 


RESULT AND DISCUSSION 


This Chapter is divided into two parts: The first part consists of validation of 
computer code, while the second part consists of numerical results for different cases. 


4,1 VALIDATION OF THE PRESENT METHOD 


In order to validate the FEM code, a static problem of a plate with a crack in the 
centre subjected to uniform tension at edges parallel to the crack axis, is solved as 
shown in Fig. 4.1a. For the analysis only one fourth of the plate is considered (due to 
the symmetry) and the mesh used is shown in Fig. 4.1b. 

The material properties and geometry are as follows ; 


Modulus of Elasticity 

E = 210 GPa 

Poissions ratio 

u =0.3 

Density 

p = 7840 Kg/m 

Length of the Plate 

L = 0.104 m 

Width of the Plate 

W = 0.04 m 

Thickness of the Plate 

B = 0.024 m 

Crack Length 

2a = 0.024 m 

Virtual crack extension 

Aa = 0.001 m 

Far field stress 

cr^ = 5 MPa 


2 

Theoretical value of energy release rate 4.084 J/m 

2 

Energy release rate using present analysis 3.967 J/m 
Thus giving percentage in error =2.8% 

Thus, it can be observed that energy release rate obtained in the present analysis is good 
agreement with theoretical result (within 3%). 
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4.2 CRACK PROPAGATION IN QUASISTATIC CASE AND DETERMINATION 
OF C: 

For this case a wedge loaded double cantilever beam specimen (DCB) is 
considered as shown in Fig. 4. 2a. Due to the symmetry, only the upper half of the 
specimen is modeled by finite elements and the mesh used is shown in Fig. 4.2b. 

The material properties and geometry are as follows: 


Modulus of Elasticity 

E = 210 Gpa 

Poissions ratio 

V = 0.3 

Density 

p = 7840 kg/m' 

Length of specimen 

L = 0.06 m 

Height of specimen 

h = 0.027 m 

Thickness of specimen 

B = 0.024 m 

Crack length 

a = 0.04 m 

prescribed displacement 

u = 0.025 mm 


As explained in the sec. 3.2, the value of stiffness Kq for the same configuration 
(geometry, bondary conditions, etc) has been obtained before the application of stiffness 
release model. 


Fig. (4.2c) shows energy release rate due to the crack propagation with increase 
in crack length, for two different values of C. It is observed that using C = 0 energy 
release rate is constant within the element but no continuous from element to element. 
After a few trials using a value of -0.045 for C it is found that energy release rate is 
linearly decreasing within the element and continuous from element to element. Thus 
it can be seen that C value can be adjusted to the correct value easily. 
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To validate the energy release rate evaluation from the energy removed in the 
decrease of stiffness release J-integral value and energy release rate value are detemined 
for the same problem. Fig. 4. 2d shows the comparison and it can be seen they match 
very well. From now onwards energy release rate G is used as it is easy to evaluate and 
far more general. A comparison of J-integral between force release model and stiffness 
release model as the crack propagates is shown in Fig. 4. 2e with increase in crack length. 

It can be seen that J-integral obtained using stiffness release model is linearly decreasing 
within the element and is continuous from element to element and J-integral obtained 
using force release model is not linearly decreasing within the element. Thus, it can be 
seen that the stiffness release model gives gives far better results than force release 
model. 

4.3 CRACK PROPAGATION AT CONSTANT SPEED IN DYNAMIC CASE 

For this case a double cantilever beam (DCB) problem considered by Verma 
(1995) is analysed. Due to the symmetry, only upper half of the specimen is modeled 
by finite elements and mesh used is shown in Fig. 4.3. Verma (1995) experimentally 
measured crack length and velocity history which were used as the input data to calculate 
the energy release rate with respect to time. Material and geometry is same as 
mentioned in the previous section. 

Fig. 4.4 presents the comparison of energy release rate due to crack propagation 
obtained using the present analysis and J-lntegral calculated by Verma (1995) using force 
release model with respect to time. The crack propagartion speed measured was of the 
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value 1750 m/sec. which is about 0.6 C„ where C„ is shear wave speed. 

3 3 


The Fig. 4.4 shows the J-integral value as a function of time evaluated by Verma 
(1995) using force release model. It may also be observed the J value was zero until t 
= 25 fxsec. Which was the time taken for the stress pulse to reach the crack tip. Fig. 
4.4 also shows the energy release rate G calculated by present method. It can be 
observed that G evaluated once the crack starts propagating. It can be seen that J- 
integral value obtained by force release model becomes highly oscillatory and unstable 
after the crack propagates by few elements. On the other hand the present method gives 
smooth values without any oscillation. Another point worthnoting is force release model 
predicts J-integral value becomes zero after opening few elements. Whereas stiffness 
release model predicts a small value while the crack is propagating at high speed. Thus 
it may be infered present model predicts more accurate and consistent is-the predictions. 

Fig. 4.5, 4.6 and 4.7 shows the similar studies make on other experimental data. 
Comparision of G as a function of crack propagation speed can be seen in the following 
table. 


Expt. No 

Propagation Speed 

Propagation Toughness 

m/s 

J/m^ 

1 

1750 

41 

2 

1500 

26 

3 

1350 

40 

4 

1200 

39 
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It appears the propagation toughness increases, but very slowly with crack 
propagation speed. Any stronger conclusion requires more data and more accurate 
measurements. 


4.4 DESIGN PROBLEM - PREDICTION OF CRACK PROPAGATION: 

I’he present method can be easily applied to predict the crack propagation history 
under the known loads, if crack propagation resistance of the material is known. That 
means that fracture resistance of the material, G^. as a function crack propagation 
velocity is supposed to be a known input data. Under a hypothetical data a limited study 
has been made, the geometry and loading are same as section 4.3. Under the 
assumption of crack initiation at time t = 40 ixsec. further crack history is evaluated. 
Fig. 4.8 shows time vs.’extension of crack and Fig. 4.9 shows propagation speed as a 
function of time. It should be noted that if the material data is fully available the 
initiation time also can be predicted. 
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Fig. 4.1b Mesh for central crack problem (one fourth) 
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Fig. 4.2b Mesh(coarse) for DCB Specimen (one half) 
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Fig. 4.2c Variation of G with increase in crack length 
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Fig. 4.2d Comparison of G and J with increase in crack length 
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Fig. 4.2e Coniparision of J between present model and force release model 
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Fig. 4.4 Comparison of G and J with time for Expt.l 
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Energy release rate, J— integral (J/ 



Fig. 4.5 Comparison of G and J with time for Expt. 2 
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Fig. 4.6 Comparison of G and J with time for Expt. 3 
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Fig. 4.8 Variation of crack extension with time 
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CHAPTERS 


CONCLUSIONS AND SCOPE FOR FUTURE WORk 

5.1 CONCLUSIONS: 

The present work proposed a stiffness release model to investigate high speed 
crack propagation under dynamic conditions. Based on the results and discussions 
present in Chapter4 the following conclusion are drawn : 

1. Stiffness release model offers non oscillatory and more stable solution. 

2. It predicts a finite (though small) propagation toughness unlike force 
release model. 

3. Design problems particularly those involving caack arrest can be easily 
solved. 

5.2 SCOPE FOR FUl URE WORK: 

1 . Any strong conclusions on propagation toughness require more 
experimental data and accurate measurements. 

2. The predetermination of Kq value by quasitatic analysis may require some 
study so as to simplfy procedure. 
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